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Abstract 

We introduce in this document a direct method allowing to solve numerically inverse type 
problems for linear hyperbolic equations. We first consider the reconstruction of the full 
solution of the wave equation posed in II x (0, T) - Q a bounded subset of lA - from a partial 
distributed observation. We employ a least-squares technique and minimize the lAnorm of 
the distance from the observation to any solution. Taking the hyperbolic equation as the 
main constraint of the problem, the optimality conditions are reduced to a mixed formulation 
involving both the state to reconstruct and a Lagrange multiplier. Under usual geometric 
optic conditions, we show the well-posedness of this mixed formulation (in particular the 
inf-sup condition) and then introduce a numerical approximation based on space-time finite 
elements discretization. We prove the strong convergence of the approximation and then 
discussed several examples for N = 1 and N = 2. The problem of the reconstruction of both 
the state and the source term is also addressed. 
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tion 
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1 Introduction - Inverse problems for hyperbolic equations 

Let LI be a bounded domain of M. N ( N > 1) whose boundary <917 is Lipschitz and let T > 0. We 
note Qt := II x (0,T) and E t := <917 x (0,T). We are concerned in this work with inverse type 
problems for linear hyperbolic equation of the following type 

f ytt-'V-(c(x)Vy) + d(x,t)y = f, (x, t) G Qt 

< y = o, (x,t)er T (i) 

l (y(-,0),y t (-,0)) = (y 0 ,yi), xeCi. 

We assume that c G C ,1 (II,K) with c( x) > cq > 0 in II, d e L°°(Qt), (yo,yi) & H := L 2 (II) x 
H- x ( II) and f &X := L 2 (0, T; ff -1 (I7)). 

For any (yo,yi) € H and any / G A', there exists exactly one solution y to ([!]), with y G 
C ,o ([0,r];L 2 (H))nC 1 ([0,T];ff" 1 (II)) (see [ID]). 

In the sequel, for simplicity, we shall use the following notation: 


Ly :=j/«-V • (c(x)Vy)+d(x,t)y. (2) 

*Laboratoire de Mathematiques, Universite Blaise Pascal (Clermont-Ferrand 2), UMR CNRS 6620, Campus de 
Cezeaux, 63177, Aubiere, France. E-mails: nicolae.cindea@math.univ-bpclermont.fr. 

'Laboratoire de Mathematiques, Universite Blaise Pascal (Clermont-Ferrand 2), UMR CNRS 6620, Campus de 
Cezeaux, 63177, Aubiere, France. E-mail: arnaud.munch@math.univ-bpclermont.fr. 


1 



1 INTRODUCTION - INVERSE PROBLEMS FOR HYPERBOLIC EQUATIONS 


2 


and X' := L 2 (0,T; H£(Q)). 

Let now w be any non empty open subset of Ll and let qx := w x (0, T) C Qx- A typical inverse 
problem for (Jlj) is the following one : from an observation or measurement y a b s in L 2 (qx) on the 
sub-domain qx, we want to recover a solution y of the boundary value problem 0 which coincides 
with the observation on qx- 

Introducing the operator P : L 2 (Qx) — > X x L 2 (qx) defined by Py := (Ly,y\ qT ), the problem 
is reformulated as : 

find y G L 2 (Q T ) solution of P y = (f,y obs ). (IP) 

From the unique continuation property for ([!]), if the set qx satisfies some geometric conditions 
and if y obs is a restriction to qx of a solution of 0 . then the problem is well-posed in the sense 
that the state y corresponding to the pair (y 0 bs , /) is unique. 

In view of the unavoidable uncertainties on the data y obs (coming from measurements, numerical 
approximations, etc), the problem needs to be relaxed. In this respect, the most natural (and widely 
used in practice) approach consists to introduce the following extremal problem (of least-squares 
type) 

minimize over H J(y 0 ,yi) ■= ^11 V ~ Vobs \\h {qT) 
where y solves 0, 

since y is uniquely and fully determined from / and the data (yo, yfi). Here the constraint y—y 0 b s = 


{ 


0 in L 2 (qx) is relaxed; however, if y obs is a restriction to qx of a solution of (llj), then problems (LS) 


and (IP) obviously coincide. A minimizing sequence for J in H is easily defined in term of the 


solution of an auxiliary adjoint problem. Apart from a possible low decrease of the sequence near 
extrema, the main drawback, when one wants to prove the convergence of a discrete approximation 
is that, it is in general not possible to minimize over a discrete subspace of {y; Ly— / = 0} subject 
to the equality (in X) Ly — / = 0. Therefore, the minimization procedure first requires the 
discretization of the functional J and of the system ([!]) ; this raised the issue of uniform coercivity 
property (typically here some uniform discrete observability inequality for the adjoint solution) 
of the discrete functional with respect to the approximation parameter. As far as we know, this 
delicate issue has received answers only for specific and somehow academic situations (uniform 
Cartesian approximation of LI, constant coefficients in 0). We refer to ns na in na and the 
references therein. 


More recently, a different method to solve inverse type problems like (IP) has emerged and 


use so called Luenberger type observers: this consists in defining, from the observation on qx, 
an auxiliary boundary value problem whose solution possesses the same asymptotic behavior in 
time than the solution of ([I]): the use of the reversibility of the hyperbolic equation then allows 
to reconstruct the initial data (yo,yi)- We refer to [8J 23] and the references therein. But, for 
the same reasons, on a numerically point of view, these method require to prove uniform discrete 
observability properties. 

In a series of works, Klibanov and co-workers use different approaches to solve inverse problems 
(we refer to [T5] and the references therein): they advocate in particular the quasi-reversibility 
method which reads as follows : for any e > 0, find y E G A the solution of 


(PVet Py) XxL 2 (q T ) + e (Vet y) A ~ ((ft Vobs)t Py) X' xL 2 (q T ),X xL 2 (q T ) ’ 


(QR) 


for all y G A, where A denotes a Hilbert space subset of L 2 (Qx) so that Py G X x L 2 (qx) for 
all y G A and e > 0 a Tikhonov like parameter which ensures the well-posedness. We refer for 
instance to Il3l where the lateral Cauchy problem for the wave equation with non constant diffusion 


is addressed within this method. Remark that (QR ) can be viewed as a least-squares problem since 
the solution y e minimizes over A the functional y —► ||Py — (/, y ob s) llxxi 2 (q T ) Eventually, 

if yobs is a restriction to qx of a solution of dTl), the corresponding y e converges in L 2 (Qx) toward 
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to the solution of (IP) as e —> 0. There, unlike in Problem (LS), the unknown is the state variable 
y itself (as it is natural for elliptic equations) so that any standard numerical methods based on a 
conformal approximation of the space A together with appropriate observability inequalities allow 
to obtain a convergent approximation of the solution. In particular, there is no need to prove 
discrete observability inequalities. We refer to the book [2]. We also mention mm where a similar 
technique has been used recently to solve the inverse obstacle problem associated to the Laplace 
equation, which consists in finding an interior obstacle from boundary Cauchy data. 

In the spirit of the works mmm, we explore the direct resolution of the optimality conditions 
associated to the extremal problem (LS), without Tikhonov parameter while keeping y as the 
unknown of the problem. This strategy, which avoids any iterative process, has been successfully 
applied in the closed context of the exact controllability of ([l]) in [TJj and mm- The idea is to 
take into account the state constraint Ly — f = 0 with a Lagrange multiplier. This allows to derive 
explicitly the optimality systems associated to (LS I in term of an elliptic mixed formulation and 
therefore reformulate the original problem. Well-posedness of such new formulation is related to 
an observability inequality for the homogeneous solution of the hyperbolic equation. 

The outline of this paper is as follow. In Section [2] we consider the least-squares problem © 
and reconstruct the solution of the wave equation from a partial observation localized on a subset 
qt of Qt ■ For that, in Section 2.1 we associate to © the equivalent mixed formulation 0 which 
relies on the optimality conditions of the problem. Assuming that qt satisfies the classical geometric 
optic condition (Hypothesis 1, see ©» , we then show the well-posedness of this mixed formulation, 
in particular, we check the Babuska-Brezzi inf-sup condition (see Theorem 2.1). Interestingly, in 


Section [2T2| we also derive a equivalent dual extremal problem, which reduces the determination 
of the state y to the minimization of an elliptic functional with respect to the Lagrange multiplier. 
In Section [3j we apply the same procedure to recover from a partial observation both the state 
and the source term. Section [4] is devoted to the numerical approximation, through a conformal 
space-time finite element discretization. The strong convergence of the approximation (yh, fh) is 
shown as the discretization parameter h tends to zero. In particular, we discuss the discrete inf-sup 
property of the mixed formulation. We present numerical experiments in Section [5] for LI = (0,1) 
and 0 C R 2 , in agreement with the theoretical part. We consider in particular time dependent 
observation zones. Section [6] concludes with some perspectives. 


2 Recovering the solution from a partial observation: a 
mixed re-formulation of the problem 


In this section, assuming that the initial (yo, yi) G H are unknown, we address the inverse problem 
(IP). Without loss of generality, in view of the linearity of the system ([I]), we assume that the 
source term / = 0. 

We consider the non empty vectorial space Z defined by 


Z:={y:y£ C([0, T], L 2 (H)) n C^QO, T], H~Ly € X}. 
and then introduce the following hypothesis : 


(3) 


Hypothesis 1 There exists a constant C 0 b s = C{lo,T, ||c|| C i^, ||d||z,=°(n)) such that the following 
estimate holds : 

||y(-,0),y t (-,0)|||j < C ohs (\\y\\ 2 L 2 {qT) + \\Ly\\ 2 x ^j, \/y £ Z. {%) 

Condition © is a generalized observability inequality for the solution of the hyperbolic equa¬ 
tion: for constant coefficients, this estimate is known to hold if the triplet (ui,T,Ll) satisfies a 
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geometric optic condition. We refer to pQ. In particular, T should be large enough. Upon the 
same condition, pTj ) also holds in the non-cylindrical situation where the domain u> varies with 
respect to the time variable: we refer to [7] for the one dimensional case. For non constant velocity 
c and potential d, we refer to m and the references therein. 

Then, within this hypothesis, for any rj > 0, we define on Z the bilinear form 


<■ y,y)z-= yydxdt + y {Ly,Ly) H -i {n) dt Vy,y G Z. 

J J qT J 0 


(4) 


In view of ©. this bilinear form defines a scalar product over Z. Moreover, endowed to this 
scalar product, we easily obtain that Z is a Hilbert space (see [7], Corollary 2.4). We note the 
corresponding norm by \\y\\ z := sj (y, y)z■ 

Then, we consider the following extremal problem : 


inf J(y) := 7j\\y ~ yoba\\ 2 L 2 ( qT ), 
subject to y £ W 


(-P) 


where W is the closed subspace of Z defined by 


W := {y G Z- Ly = 0 in X} 

and endowed with the norm of Z. 

The extremal problem © is well posed : the functional J is continuous over W, is strictly 
convex and is such that J(y) —► +oo as ||y||w —► oo. Note also that the solution of 0 in W does 
not depend on y. 

Remind that from the definition of Z , Ly belongs to X. Similarly, the uniqueness of the solution 
is lost if the hypothesis © is not fulfilled, for instance if T is not large enough. Eventually, from 
( |W| ), the solution y in Z of 0 satisfies (y(-, 0), y t (-, 0)) € H, so that problem 0 is equivalent to 
the minimization of J with respect to (ycnl/i) G H as in problem (IP), Section 1. 


We also recall that for any z G Z there exists a positive constant Ciyr such that 


(5) 


( 6 ) 


IMIluqt) < 0), z t (-, OjUjj + ||F'^||x 

This equality and © imply that 

INI L 2 (Qt) — Ctl,T ^Cobs||~||| 2 ( gT ) + (1 + C 0 bs)\\Lz\\ 2 ^ , \/z G Z. 

2.1 Direct approach 

In order to solve 0 , we have to deal with the constraint equality which appears in the space 
W. Proceeding as in [12], we introduce a Lagrangian multiplier A G X' and the following mixed 
formulation: find (y, A) G Z x X' solution of 


a(y, y) + &(y, A) = l(y), Vy g Z 

b(y,X) = 0, VAGI', 


where 


a : Z x Z — > R, a(y,y):= // yydxdt, 

J J qT 

b-.ZxX'—tM., b(y,X):= f (A, Ly) H u n ^ H -i^dt, 

Jo 

l ■■ Z M, l(y) := y obs y dxdt. 

J J qx 


(7) 

( 8 ) 
(9) 

( 10 ) 
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System ([7| is nothing else than the optimality system corresponding to the extremal problem Im¬ 
precisely, the following result holds : 

Theorem 2.1 Under the hypothesis 

(i) The mixed formulation 0 is well-posed. 

(ii) The unique solution (y, A) £ Z x X' to ([7|| is the unique saddle-point of the Lagrangian 
C : Z x X' —► M defined by 

£{y, A) -=^a{y, y) + b(y, A) - l(y). 

(Hi) We have the estimate 

\\y\\z = \\y\\L 2 (q T ) < \\yobs\\L 2 (q T ), I|A||.y' < + v\\yob S \\L^( qT )- (u) 


PROOF- We use classical results for saddle point problems (see [4], chapter 4). 

We easily obtain the continuity of the bilinear form a over Z x Z, the continuity of bilinear b 
over Z x X' and the continuity of the linear form l over Z. In particular, we get 

\\l\\z> = \\yobs\\L 2 (q T ), IMI (zxzy = 1, \\b\\ (zxx'y = y~ 1/2 - (12) 

Moreover, the kernel M(b) = {y £ Z\ b(y, A) = 0 VA G X'} coincides with IV: we easily get 

a (y,y) = h\\ 2 z, Vy e W(6) = w. 


Therefore, in view of [H Theorem 4.2.2], it remains to check the inf-sup constant property : 3<5 > 0 
such that 

b{y, A) 


inf sup 


> < 5 . 


\ex’ y&z ||y||z||A||x' 

We proceed as follows. For any fixed A £ X', we define y as the unique solution of 
Ly = -AA in Q Tl (j/(-,0),y t (-,0)) = (0,0) on fi, y = 0 on S T . 
We get b(y, A) = ||A|| 2 V , and 

Iblll = \\y\\h( qT )+y\\M\x'- 

Using ([ 5 ]), the estimate \\y\\L 2 (q T ) < x/Un.rllA||x' implies that y £ Z and that 

b(y, A) _ 1 


sup 


> 


yez ||y||z||A|U' - VCn,T + y 


> 0 


leading to the result with 5 = (Cn,T + v) 1 ^ 2 - 

The third point is the consequence of classical estimates (see [4], Theorem 4.2.3.) 


\z < —\\l\\z', 
OL 0 


IX' < \ ( 1 + ^ 
«o 


where 


ao := inf 


5 

a(y,y) 


yew) inn ' 


(13) 


(14) 


(15) 


Estimates (12) and the equality ao = 1 lead to the results. Eventually, from (12), we obtain that 

2 

l|A||.Y' < ^||2/o&s|U 2 (g T ) 
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and that i5 > (Cq^t + v) to get (111. □ 

In practice, it is very convenient to ’’augment” the Lagrangian (see m) and consider instead 
the Lagrangian C r defined for any r > 0 by 


£ r {y, A) := ^ a r (y , y) + b(y, A) - /(y), 
a r (y,y) ■■= a(y,y) +r\\Ly\\ 2 x . 


Since a r (y,y ) = a(y,y) on W, the Lagrangian C and C r share the same saddle-point. The positive 
number r is an augmentation parameter. 


Remark 1 Assuming additional hypotheses on the regularity of the solution A, precisely LX £ 
L 2 (Qt) and (A,At)|t=o,r £ Hq(H) x L 2 (U), we easily prove, writing the optimality condition for 
£, that the multiplier A satisfies the following relations : 

f LX = -(y - y obs ) l u in Q T , A = 0 in E T , 

< 16) 
[X = Xt = 0 on O x {0, T}. 


Therefore, X (defined in the weak sense) is an exact controlled solution of the wave equation through 
the control -{y - y a bs ) 1^ £ L 2 (q T ). 

• If dobs is the restriction to qr of a solution of then the unique multiplier X must van¬ 
ish almost everywhere. In that case, we have sup AgA inf^gy C r (y, A) = inf ye v £ r (y, 0) = 
inf ye v Jriy) with 

Jr(y) '■= 2 II y - 2/o&s||l 2(Q t ) + ^ll^yllx- ( 17 ) 

The corresponding variational formulation is then : find y £ Z such that 


a r {y,y) 


// yydxdt + r (Ly, Ly) H -i( Q) dt = l(y), 
J J Qrp J 0 


Vy £ Z. 


• In the general case, the mixed formulation can be rewritten as follows: find ( z , A) £ Z x X' 
solution of 


( Pry,P r y}xxL 2 (q T ) + {Ly, X)x,X' = ((0 ,yobs), Pry}xxL 2 (q T ), Vy £ Z, 
{Ly , X)x,x' =0, VA £ X' 


(18) 


with P r y := (y/rLy,y\ qT ). This approach may be seen as generalization of the (QR) problem 
(see (QR)), where the variable X is adjusted automatically (while the choice of the parameter 


e in (QR) is in general a delicate issue). 


System (16) can be used to define a equivalent saddle-point formulation, very suitable at the 
numerical level. Precisely, we introduce - in view of (16) - the space A by 


A := {A : A £ C([0, TfH^n)) n ^([0, T]; L 2 (L>)), 

LX £ L 2 (Qt), X(-,0) = A t (-,0) = 0}. 


Endowed with the scalar product (A, A)a := 


(A A + LXLX) dxdt, we check that A is a Hilbert 


' Q-i 


space. Then, for any parameter a £ (0,1), we consider the following mixed formulation : find 
(y, A) £ Z x A such that 
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(19) 


where 


I a r , a {y, y) + b a (y, A) = h, a (y), Vy £ Z 

\ b a (y, A) - c a (X, A) = h,a(X), VX £ A, 

a^:Zx2->I, a r , a (y ,y) := {1 - a) yy dxdt + r (Ly, Ly) H -i^dt, 

J J q^p J 0 

b a :Zx A-fM, b a {y, A) := / (X, Ly) H i {a) H -i^dt - a yLXdxdt, 

J 0 J J q'j' 

c a : A x A —> R, c a (A, A) := a // LX LX, dxdt 

J J q t 

h, a :Z^R, h, a (y) := (1 - a) // y obs ydxdt, 

J J qrjp 

h,a ■ A ->• R, / 2 , a (A ) :=-a y obs LXdxdt. 

J J qrp 

From the symmetry of a rjC( and c a , we easily check that this formulation corresponds to the 
saddle point problem : 

sup inf jC rta (y, A), 

AeA vcz 

£-r,a(y> A) := E r (y, A) — — ||TA + (y — yobs)^ui Hl 2 (q t )' 


Proposition 2.1 Under the hypothesis ©. for any a £ (0,1), the formulation |Vff| ) is well-posed. 
Moreover, the unique pair (y, X) in Z x A satisfies 


0i\\y\\z + » 2 ||A||i < ^ I- l|j/o6s|||^( 9T )- 


( 20 ) 


with 


0\ \= mini 1 — a, — 
V 


92 :=1 2 min ( a 'Ch 


PROOF- We easily get the continuity of the bilinear forms a r , a , b a and c a : 

|ar,a(j/,y)| < max(l - a, -)||y||z||y||z, Vy,y £ Z, 

V 

\b a (y, A)| < max(a, -^=)||y||z||A|| A , Vy e Z,MX £ A, 

Vv 

|c a (A,A)<a||A|| A ||A|| A , VA, A e A 

and of the linear form h and l 2 ■ \\h\\z' = (1 - ot)\\y obs \\ L i( qT ) and ||Z 2 ||a' = ot\\y obs \\ L 2 ( qT ). 
Moreover, since a £ (0,1), we also obtain the coercivity of a r , a and of c a : precisely, 

a r , a {y,y) > min ^1 - a, ^ \\y\\ 2 z , Wy £ Z, 

c a (A, A) > min( am, TO ) ||A||a VA £ A, Vm6 (0,1). 

\ Of 2 ,t J 


The result [3] Prop 4.3.1] implies the well-posedness and the estimate (20) taking m = 1/2. □ 

The a-term in £ r a is a stabilization term: it ensures a coercivity property of C r ^ a with respect 
to the variable A and automatically the well-posedness. In particular, there is no need to prove 
any inf-sup property for the application b a . 
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Proposition 2.2 If the solution (y, A) e Z x X' of 0 enjoys the property X £ A, then the 
solutions of 0 and [Hty coincide. 


PROOF- The hypothesis of regularity and the relation (16) imply that the solution (y, X) £ Z x X' 
of 0 is also a solution of (19). The result then follows from the uniqueness of the two formulations. 
□ 


2.2 Dual formulation of the extremal problem ([7]) 

As discussed at length in [T2], we may also associate to the extremal problem 0 an equivalent 
problem involving only the variable A. Again, this is particularly interesting at the numerical level. 
This requires a strictly positive augmentation parameter r. 

For any r > 0, let us define the linear operator V r from X' into X' by 

V r X := -A ~ 1 (Ly), VA e X' 

where y £ Z is the unique solution to 

a r (y,y) = b{y, A), My £ Z. (21) 


The assumption r > 0 is necessary here in order to guarantee the well-posedness of (21). Precisely, 
for any r > 0, the form a r defines a norm equivalent to the norm on Z. 

The following important lemma holds: 

Lemma 2.1 For any r > 0, the operator V r is a strongly elliptic, symmetric isomorphism from X' 
into X'. 

PROOF- From the definition of a r , we easily get that ||7VA||.x' < r —1 1| A||jsf/ and the continuity of 


V r - Next, consider any A' € X' and denote by y' the corresponding unique solution of (21) so that 
V r X' := — A ~ 1 (Ly'). Relation (21) with y = y' then implies that 


/ ('PrX l ,X) H u n) dt = a r (y,y l ) 
J 0 


( 22 ) 


and therefore the symmetry and positivity of V r . The last relation with X' = X and the observability 
estimate © imply that V r is also positive definite. 

Finally, let us check the strong ellipticity of V r , equivalently that the bilinear functional 
(A, A') —> fo(V r A, A , )#i(Q) iff i(Q) dt is Af'-elliptic. Thus we want to show that 


/VrA, A>*i (n) dt > C||A|| 2 y „ VA € X' 

J 0 


(23) 


for some positive constant C. Suppose that (23) does not hold; there exists then a sequence 
{A„} n >o of X' such that 


IIA 


n X‘ 


= 1, Mn > 0, lim f (V r X n , X n ) H i (n) dt = 0. 

n—too Jq u v 


Let us denote by y n the solution of (21) corresponding to A„. From (22), we then obtain that 

lim r\\Ly n f x + ||y„|| 2 2 ( } = 0. 

n- 4-oo \ tJ- J 


(24) 


From (21) with y = y n and A = A n , we have 

rT 


[ (r(-A x )Ly n - A„, (-A 1 )Ly) Rl dt+ If y n ydxdt = 0, My £ Z. 
J 0 0 J J qt 


( 25 ) 
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We define the sequence {y n } n > o as follows : 


i Ly n =rLy n + A _1 A„, in Q T , 

Vn = 0, in s t, 

¥n(->°) =y n ,t(->°) = 0> in 


so that, for all n, y n is the solution of the wave equation with zero initial data and source term 
rLy n + AA„ in X. Using again JHll, we get \\y n \\L 2 (q T ) < \fCnpr\\rLy n + AA n ||x, so that y n G Z. 


Then, using (25) with y = y n we get 

lk(—A ~ 1 )Ly n - A n ||x' < VC^\\yn\\m qT ). 


Then, from (24 1, we conclude that limn^+oo l|An||x' = 0 leading to a contradiction and to the 
strong ellipticity of the operator V r . □ 

The introduction of the operator V r is motivated by the following proposition : 


Proposition 2.3 For any r > 0, let yo G Z be the unique solution of 

a r {yo,y) = Ky), Vy&z 

and, let J** : X' —> X' be the functional defined by 

J**(A) = - J (V r \,X) H i(ii)dt — b(yo,X)- 

The following equality holds : 

sup inf C r (y, A) = - inf J**(A) + £ r (y 0 ,0). 

\GX'V^ Z x & x 

The proof is classical and we refer for instance to 121 in a similar context. This proposition 
reduces the search of y, solution of problem <0 to the minimization of J**. The well-posedness 
is a consequence of the ellipticity of the operator V r . 


Remark 2 The results of this section apply if the distributed observation on qx is replaced by a 
Neumann boundary observation on a sufficiently large subset Ey ofdflx (0,T) (i.e. assuming ^ = 
yobs G L 2 (Tit) is known on Tit)- This is due to the following generalized observability inequality: 
there exists a positive constant C 0 b s = C(w, T, ||c|| C i^, ||d||z,= (n)) such that the following estimate 
holds : 


\\y(‘i o), 2/tO, o)||# i(Q)xi 2 (fi) — Cobs 


dv 


2 

T 2 (Et) 


I Ly \\ 2 


l 2 (Q t ) h 


Vy&Z 


(26) 


which holds if the triplet (Qx,Z,x,T) satisfies the geometric condition as before (we refer to 11 Of 
and the 7'eferences therein). Actually, it suffices to re-define the form a in ^ by a(y,y) := 
ffs T dtdt dadx and the f° rm 1 b y %) := Ue t dhVobs dadx for all y,y G Z. 


Remark 3 We emphasize that the mixed formulation 0 has a structure very closed to the one 
we get when we address - using the same approach - the null controllability of more precisely, 
the control of minimal L 2 (qx)-norm which drives to rest the initial data {yo,yi) G Hq(U) x L 2 (fl) 
is given by v = <p l qT where (ip, A) G $ x f 2 (0,T; Hq (U)) solves the mixed formulation 


a(ip,p) + b(tp,\) = l (ip) 

b(y>, A) = 0, 


\/p G d> 

VA G L 2 (0,T ; Hq(S!)), 


( 27 ) 
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where 


a : $ x <f> —> R, a(ip, ip) = // ip(x,t)ip(x,t) dx dt 

J J qT 

b : $ x L~(Q,T\Hq(LI)) > R, b(ip, A) = f (Lip, X) H ~i H idt 

Jo 

l : $ —» R, l(ip) = — 0), + / p(-,ti)y\dx. 

J 0 


with $ = (p £ L 2 (Qt), ip = 0 on S x such that Lp £ L 2 (0,T;H 1 (J1))}. We refer to [Wj. 

Remark 4 Reversing the order of priority between the constraint y — y a b s = 0 in L 2 (qT) and 
Ly — / = 0 in X, a possibility could be to minimize the functional y —> ||Ly — f\\x over y £ Z 
subject to the constraint y — y 0 b s = 0 in L 2 (qT) via the introduction of a Lagrange multiplier in 
L 2 (qT). The proof of the following inf-sup property : there exists S > 0 such that 


IL 


qT A y dxdt > 


inf sup 
a £L 2 {q T )y£z ||A||z, 2 (g T )||y| 


associated to the corresponding mixed-formulation is however unclear. 
{ QR), this difficulty disappears (we refer again to the book f75j/). 


If a e-term is added as in 


3 Recovering the source and the solution from a partial ob¬ 
servation: a mixed re-formulation of the problem 


Given a partial observation y 0 b s of the solution on the subset qx C Qt, we now consider the 
reconstruction of the full solution as well as the source term / assumed in X. We assume that the 
initial data (yo,Vi) € H are unknown. 

The situation is different with respect to the previous section, since without additional assump¬ 
tion on /, the couple (y, f) is not unique. Consider the case of a source / supported in a set which 
is near dfl x (0,T) and disjoint from qT'. from the finite propagation of the solution, the source 
/ will not affect the solution y in qx- On the other hand, the determination of a couple (y, f ) 
which solves 0 such that y coincides with y obs is straightforward : it suffices to ’’extend” y on 
Qt \ qr appropriately to preserve the boundary conditions, then compute Ly and recover a source 
term. However, we emphasize that, on a practical viewpoint, the extension of y 0 bs out of qx is not 
obvious. Moreover, this strategy does not offer any control on the object /. 

We briefly show that we can apply the method developed in Section [2] which allows a robust 
reconstruction and then consider the case of uniqueness via additional condition on /. 

We assume again that © holds. We note Y := Z x X and define on Y the bilinear form, for 
any e, r] > 0 


((y,f)i(yJ))Y--=H yydxdt + rj I (Ly - f,Ly - f) H -i^dt 

/ qrp J 0 

rT 


(28) 


+ £ / (/> f)H- 1 <n)dt, v(y, /),(?/,/) G Y. 


/o 


In view of ©. this bilinear form defines a scalar product over Y. Moreover, endowed to this scalar 
product, we easily obtain that Y is a Hilbert space (we refer to ||7j). We note the corresponding 
norm by ||(y,/)||y := y/((y, /), (y, f))v- 
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Then, for any e > 0, we consider the following extremal problem : 


(Ve) 


inf J £ (y, f) := * ||y - y obs f L 2 {qT) + |||/|| 2 Y 
subject to (y, /) £ W 


where W is the closed subspace of Y defined by IT := {(y, /) £ Y ; Ly — f = 0 in X} and endowed 
with the norm of Y: precisely, it follows that 

\\(y,f)\\w ■= JWvWhM+'Wm, v(y,/)eir. 

The extremal problem ( V e ) is well posed : the functional J e is continuous over W, is strictly 
convex and is such that J £ (y, /) —► +oo as ||(y, f)\\w —> oo. Note also that the solution of (V e ) in 
W, depends on e but not on rj. 

Remark also that if e = 0, then J e is a priori only convex leading possibly to distinct minima. 
This justifies the introduction of the e-term in the functional J e . We emphasize however that the 
e-term is not a regularization term as it does not improve the regularity of the state y. 

Eventually, from (|W]), the solution (y £ , f e ) in W of (V e ) satisfies (y e (-, 0), y £j t(-, 0)) £ H, so that 
problem (V e ) is again equivalent to the minimisation of J £ with respect to (yo,yi,/) € H x X. 

Proceeding as in Section [5J we introduce a Lagrangian multiplier A £ G X' and the following 
mixed formulation: find ((y e ,/ £ ), A e ) G Y x X' solution of 


where 


ae((yeje),{yj)) + b((y,f),\ £ ) = %,/), V(y,/)Gb 

K(y e ,/ e ), a) = o, va g x', 


a E : 7 x y G I, a £ ((y,/), (y,/)) := // yy dxdt + e(f, f) x , 

J J qj- 

r T 

b :Y x X 1 

/ o 


b{{y, /), A) / (A, Ly - f) H un), H -i(n) dt, 
Jo 


l :Y ->1, l(y, f) := // y 0 b s ydxdt. 

I q T 


(29) 


(30) 

(31) 

(32) 


Theorem 3.1 Under the hypothesis ©, the following hold : 

(i) The mixed formulation f£ff| ) is well-posed. 

(ii) The unique solution ((y £ , / £ ), A £ ) G YxX' is the saddle-point of the Lagrangian L e : YxX' —> 
R defined by 

£e{{y, /), A) := *a e ((y, /), (y, /)) + 6((y, /), A) - l(y, /). 

Moreover, the pair (y £ ,/ £ ) solves the extremal problem (V e ). 

(Hi) The following estimates hold : 

Kv s Js)\\y = (ll 2 /slli 2 ( 9T )+e||/ £ || 2 ) 1/2 < Wyotsl\mq T ) (33) 

and 

||A e ||l2(q t ) < 2 a/ Cq.,t + y||y 0 6 S ||i2( gT ) 

for some constant Cq.t > 0. 


(34) 
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The proof is very closed to the proof of Theorem 2.1 In particular, the obtention of the inf-sup 

tant 

(35) 


property is obtained by taking, for any A £ A'', / = 0 and y as in (14) so that the inf-sup constant 

b{(y,f),x) 


Se '- 

is bounded by above by (Cn,T + y) _1//2 uniformly with respect to e. 

Remark in particular that the inequality (33) implies that, at the optimality, since e > 0, the 
equality ||y - y 0 bs\\L^{q T ) = 0 can not hold if f e ^ 0. 

Remark 5 We may also prove the inf-sup property using the variable f: for any A £ X', we set 
y = 0 and f = A A £ X. We get 


sup 


K(yJ),x) 


> 


b(( 0,AA),A) 


1 


{y,f)£Y ||(y I /)IM|A||x' 11(0, AA)||v||A||x' V £ + v 

so that 5 e > (e + ? 7 ) -1 / 2 . Therefore, the estimate 


I-Mix' < -H|2 /o6 S ||l 2 (,? t ) 


implies that 


ll^ellx' < 2V £ + 'n\\yobs\\L 2 (q T )- 

This argument is valid if and only if f is distributed everywhere in Qt- 


(36) 

□ 


Remark 6 The estimate (36) implies that the multiplier X e vanishes in X' as e + y — > 0 + (recall 
that e and y can be chosen arbitrarily small in W- 

Remark 7 

(a) Assuming enough regularity on the solution X e , precisely that LX e £ L 2 {Qt) and (A, X t )t=o,T £ 
Hq(LI) x L 2 (LL), we easily check that the multiplier X e satisfies the following relations : 

LX e (y £ yobs ), Uy € f £ = 0, cf e -t- AA £ =0 in Qt> 

X e = 0 in Tit , 

A £ = A £i t = 0 on LI x {0,T}. 

Therefore, A £ is an exact controlled solution of the wave equation through the control — (y £ — 


hobs) lu and from (36) implies that 

he - y a bs\\L2(q T ) ->• 0 as £ —» 0 + . 


(37) 


Remark however that f e may not be bounded in X' uniformly w.r.t. £ (contrarily to the 
sequence {y/ef e ) E >o)- 

(b) The equality Ly e = f e becomes eLy e = — AA E and leads to L(eA~ 1 Ly e ) = —LX e = ( y e — 
hobs) i„- Finally, y e solves, at least in V, the boundary value problem 

T(£"( A )Ly e ) T y e 1^ — yobs 1 uj , in Qt , 

{sLy e ) = (£Ly e ) t = 0, in U x {0,T} 
y e = 0, on T t 

or equivalently the variational formulation: find y e £ Z (see solution of 


■ / (Ly E ,Ly) H -i( a) dt+ // y e ydxdt= // y obs ydxdt , Vy 
J 0 J J Qt «/ J Qt 


£ Z 


(38) 
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which actually can be obtained directly from the cost J E , replacing from the beginning f by 
the term Ly. From the Lax-Milgram lemma, (38) is well-posed and the following estimates 
hold : 

\\ye\\L 2 ( qT ) < \\yobs\\L 2 ( qT ), Ve\\Ly s \\x < \\yobs\\L*( qT )- 


This kind of variational formulation involving the fourth order term Ly e Ly has been derived 
and used in m in a controllability context. 


For any e > 0 and any y 0 j, s £ L 2 (qT), the method allows to recover a couple (y E , f e ) such that 
Ly e = f e in Qt and y e is closed to y 0 b s (see ([37])) . In view of the loss of uniqueness, we have no 
information on the limit of the sequence as e —> 0: the sequence may be unbounded at the limit in 
L 2 (Qt) x L 2 (Qt) even if y a b s is the restriction to qr of a solution of |T]). 

Remark 8 Contrarily to the inf-sup property, the coercivity of a e over N(b) does not hold uni¬ 
formly with respect to e. Recall that the e-term has been introduced to get a norm for Y. This 
enforces us to add this term in the mixed formulation. 


Remark 9 A fortiori, if the initial condition (yo,yi) £ H is known, one may recover the pair 
(y,f) £ Y from y 0 bs and (yo,yi). The procedure is similar; it suffices to define two additional 
Lagrange multipliers (Ai, A 2 ) £ L 2 (I2) x Hq(D) to deal with the constraint y(-, 0) = yo and yt{-, 0) = 
2 /i respectively. The extremal problem is now : 

, urf J e(y,f) ■= \h ~ yobs\\h (qT ) + ^\\f\\x> 

(v,f)€W Z Z 

where W is the closed subspace ofY defined by 

W := {( y,f ) £ Y; Ly - f = 0 in X',(y(-,0),y t (-,0)) = (y 0 ,yi) in H}. 


The corresponding mixed formulation is : find {{y e , fe), (A e , A £j i, A £i 2)) £ Y x A solution of 


r((ye,/ £ ),( 2 /,/)) + fe(( 2 /,/),(A £ ,A ei i,A £i2 )) = h(yj), 


v(y,f)&Y 


b{(y s , f s ), (A, Ai, A 2 )) — ( 2 (A, Ai,A 2 ), V(A, Ai, A 2 ) e A, 


(39) 


where a E is given by (301 and 

b:Y x A-»R, b((y,f),(X,X 1 ,X 2 )):=[f X(Ly - f) dxdt 

J J Qx 1 

+ (y('i 0), ^1 }l 2 (Q) + {yt{-, 0), A 2 ) H -i( n ) ifl -i(Q) 
h : Y -)• R, h{y,f):= ff y 0 b s ydxdt 

J J qrp 

l 2 : A —> K, l 2 { A, Ai, A 2 ) := (y 0 , Xi) L 2 ^ + (2/1, A 2 )n-i(n),_f/i(n) 

with A := X' x L 2 (LL) x Hq(Q). Using the estimate we easily show that this formulation is 
well-posed. □ 


In view of Remark [ 7 ] (a), we may also associate to the mixed formulation (29) a stabilized 
version, similarly to ©. 

Again, it is very convenient to ’’augment” the Lagrangian (see [16]) and consider instead the 
Lagrangian C e r defined for any r > 0 by 


Ar,r((y, /), A) := ^a Eir ((y, /), (2 1 , /)) + b(y , A) - l(y, /), 
a e ,ri(y, /), {y, /)) := a e {{y, /), (y, /)) + r\\Ly - f\\ 2 x . 
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Since a e (y,y) = a e ^ r (y,y) on W, the Lagrangian C e and C £ . r share the same saddle-point. The 
positive number r is an augmentation parameter. Similarly, proceeding as in Section |2.2[ we may 
also associate to the saddle-point problem sup AgX / ini/ y f )eY C re ((y, /), A) a dual problem, which 
again reduces the search of the couple (y s , f e ), solution of problem (V e ), to the minimization of a 
elliptic functional in A e . 

Proposition 3.1 For any r > 0, let ( yo, fo) GY be the unique solution of 

aeA(yo,fo),(yJ)) = v(y,7) e y 

and let V e ,r be the strongly elliptic and symmetric operator from X' into X' defined by V £<r X '■= 
—A~ 1 (Ly — f) where ( y,f) £ Y is the unique solution to 

aeA(yJ)’(yJ)) = b ({y>f)A), v(j/J)ey; (40) 

Then, the following equality holds 

sup inf £ e ,r((y, /), A) = - inf J**( A) + £ e , r (( y 0 , fo), 0). 

A EX' (!/J )0 A eA' 

where J** : X' —> X' is the functional defined by 

J** r W = 2 J (A,A, dt — b((yo, fo), A). 

Compared to the previous section, the additional unknown f e on the problem guarantees that 
the term \\y e — y 0 bs\\L 2 (q T ) vanishes at the limit in e, for any y a b s £ L 2 (qT), be a restriction of 
a solution of ([!]) or not. The situation is different if additional assumption on / enforces the 
uniqueness of the pair (y,f) (we refer to S3: and the references therein). 


4 Numerical Analysis of the mixed formulations 

4.1 Numerical approximation of the mixed formulation ([7]) 

We consider the numerical analysis of the mixed formulation ([7|, assuming r > 0. We follow US] , 
to which we refer for the details. 

Let Zh and A^ be two finite dimensional spaces parametrized by the variable h such that 
Zh C Z, Ah C X' for every h > 0. Then, we can introduce the following approximated problems: 
find the (y^, A h) £ Zh x Ah solution of 


ar{yh,y h ) + b(y h ,Xh) 
b(yh, X h ) 


KyJ, 

o, 


7 Vh e Zh 

VA h £ Ah. 


(41) 


The well-posedness of this mixed formulation is again a consequence of two properties: the coer- 
civity of the bilinear form a r on the subset 


ATh(b) = {yh £ Zh ; b(yh, Xh) = 0 VA h £ A^}. 

Actually, from the relation a r (y,y) > (r/ 77 )||J/1|^ for all y £ Z, the form a r is coercive on the full 
space Z, and so a fortiori on Afh{b) C Z\, C Z. The second property is a discrete inf-sup condition. 
We note Sh > 0 by 

r • f b(y h , X h ) 

d h := ml sup —— n — n —n—. 

\h€A h yhGZh \\Xh\\x'\\yh\\z 


( 42 ) 
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For any fixed h, the spaces Zh and A h are of finite dimension so that the infhnum and supremum 
in (42) are reached: moreover, from the property of the bilinear form a r , it is standard to check 
that Sh is strictly positive. Consequently, for any fixed h > 0, there exists a unique couple (yh, Xh) 
solution of (41). 

We then have the following estimate. 


Proposition 4.1 Let h > 0. Let (y, A) and (yh, A h) be the solution of 0 and of m respectively. 
Let Sh the discrete inf-sup constant defined by \4ty - Then, 

II y - yh\\z < 2 ^i + —0-^ d (y> z h) + 0rf(A, A/i )> ( 43 ) 

||A-A/j||x'< (^2 4— ——\—d(y,Z h )-\ — —-d(X,A h ) (44) 

V Vy°h) o h y/yOh 

where d(X,A h ) := inf Aft eA fe ||A - X h \\x> and 


d(y, Z h ) 


:= inf \\y-y h \\z 

Vh£Z h 

= y ^ Zh (llz/ - 2 ^IIl 2 (9t) + v\\L(y - yh)\\]^j 


PROOF- From the classical theory of approximation of saddle point problems (see [H Theorem 
5.2.2]) we have that 


to - „„|U < I ^ >- | d{y , Zh) 


vl|6|l 


+ 


OL o 

Hl(ZxX')' 

«0 


afidh 


d( A, A h ) 


and 


(45) 


IA — A^II.y' < 


2 ll a ^ll(ZxZ)' |K-||(ZxZ)'|HI(XxX')' 

a^S h 

3 ll a r|| 2 |H|(ZxX')' 


% 


d(y, Z h ) 


an s, 


d( A, A 


h ■ 


0 Oh 


(46) 


Since, \\a r \\(z x zy = a 0 = 1; IHI(zxA)' = the result follows. □ 

Remark 10 For r = 0, the discrete mixed formulation m is not well-posed over Zh x A^ because 
the form a r= o is not coercive over the discrete kernel ofb: the equality b(yh, A h) = 0 for all Xh £ A^, 
does not imply in general that Lyh vanishes. Therefore, the term r||T?/ft ,||which appears in the 
Lagrangian C r , may be understood as a stabilization term: for any h > 0, it ensures the uniform 
coercivity of the form a r and vanishes at the limit in h. We also emphasize that this term is not 
a regularization term as it does not add any regularity on the solution yh- 


Let nh = dim Zh,mh = dim A;, and let the real matrices A r h £ M. nh ' nh , Bh £ W 7lh ’ nh , Jh £ 
and Lh £ R n ' 1 be defined by 


ar(lJh,yh) — (A r ^h{yh} , {yh})M. n h 1 M. n h 
b(yh,Xh) = {Bh{yh},{Xh})wi™h : wi™ h 




Qt 


dx dt — 5 


. Kyh) = (Lh, {yh})R™h 


Vyh,yh £ Z hl 
V yh £ Zh, Xh £ Ah, 

VA h, Xh £ Ah, 


(47) 


V yh £ Zh, 
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where {yh} £ R" h denotes the vector associated to yh and (•, •)R"h ) R"h the usual scalar product over 


With these notations, the problem (41) reads as follows: find {yh} £ R n, ‘ and {A/,} £ R m '* 
such that 


Ar : h 

B h 


Bl 

0 


n h+ n 




{yh} 

{A4 


ri h+ rt 


L h 

0 


(48) 




The matrix A r .h as well as the mass matrix Jh are symmetric and positive definite for any h > 0 
and any r > 0. On the other hand, the matrix of order rrih + nh in (48) is symmetric but not 
positive definite. We use exact integration methods developed in m for the evaluation of the 
coefficients of the matrices. The system (48) is solved using the direct LU decomposition method. 


4.1.1 C' 1 -finite elements and order of convergence for N = 1 

The finite dimensional and conformal space Zh must be chosen such that Lyh belongs to X = 
L 2 (0,T; 7f _1 (f2)) for any yh £ Zh- This is guaranteed, for instance, as soon as iph possesses 
second-order derivatives in L 2 oc (Qt). As in [12], we consider a conformal approximation based on 
functions continuously differentiable with respect to both variables x and t. 

We introduce a triangulation 7), such that Qt = UkgT h B and we assume that {Th}h>o is a 
regular family. We note h := max{diam(7T), K £ 7*.}, where diam(lv) denotes the diameter of K. 
Then, we introduce the space Zh as follows : 

Z h = {y h £ZcC 1 (Q^):z h \ K £P(K) VK £ T h , z h = 0 on E T }, (49) 

where P(7v) denotes an appropriate space of functions in x and t. In this work, we consider two 
choices, in the one-dimensional setting (for which 12 c R, Qt CK 2 ): 

(i) The Bogner-Fox-Schmit (BFS for short) C 1 -element defined for rectangles. It involves 16 
degrees of freedom, namely the values of yh,yh,x,yh,uyh,xt on the four vertices of each rect¬ 
angle K. Therefore P [K) = P 3 iX <g> P 3 jt where P r ^ is by definition the space of polynomial 
functions of order r in the variable £. We refer to EJ ch. II, sec. 9, p. 94], 

(ii) The reduced Hsieh-Clough- Tocher (HCT for short) C 1 -element defined for triangles. This is 
a so-called composite finite element and involves 9 degrees of freedom, namely, the values of 
Vh,yh,x,yh,t on the three vertices of each triangle K. We refer to [9j ch. VII, sec. 46, p. 285] 
and to El 125 where the implementation is discussed. 

We also define the finite dimensional space 

A h = {X h £ C°(Qt), X h \K £ Q(A') VTsT £ %}■ 

where Q(K) denotes the space of affine functions both in x and t on the element K. 

For any h > 0, we have Zh C Z and Ah C X'. 

We then have the following result: 

Proposition 4.2 (BFS element for N = 1 - Rate of convergence for the 
h > 0, let k £ {1,2} be a positive integer. Let (y, A) and (yh,Xh) be the solution of 
respectively. If the solution (y, A) belongs to H k+2 {QT) x H k \Qx), then there exists 
constants 

Ki = ATj(|[y|| H k + 2 (Q T ), IMIcRQt)’ IMIU°°(Q t )), * £ {1 j 2}, 

independent of h, such that 

\\y ~ Vh\\z < K 1 ^-—^{y/r]+ ^-)(h 3 + y/rjh) + 1^ , (50) 

h k ~ 1 / 1 \ 

IIA — A/j||x' < I<2 yW 1 ! + j^)(^ 3 + \/nh) + lj . 


norm Z) Let 
{?[) and J7~/| ) 
two positives 


( 51 ) 
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PROOF - From 0 ch. Ill, sec. 17], for any A £ H k (QT), k < 2, there exists C\ = Ci(|| A||#fc(Q T )) 
such that 

||A-n Afi>r ,(A)|U' <C^h k -\ V/i > 0 (52) 

where Hk h ,T h designates the interpolant operator from X' to A h associated to the regular mesh 
Th- Similarly, for any y € H k+2 (QT ), there exist C 2 = C2(\\y\\H k + 2 (Q T )) suc ^ that for every h > 0 
we have 

\\y-H Zh ,r h {y)\\L^Q T ) < C 2 h k+2 , \\y - H Zht r h (y)\\m(Q T ) < C 2 h k . (53) 

Then, observing that 


II L V - Ly h \\ x < A'(||c|| cl( g^ } , ||d||L~ ( Q T) )||y - yh\\m(Q T ), 
for some positive constant K , we get that 

d(y, Z h ) = ^inf ^ (j|y - y h \\\i {qT) + y\\Ly - Ly h \\ 2 x ^j 

1/2 

< C 2 ( (h k+2 ) 2 +r)K 2 (h k ) 2 


(54) 


(55) 


< C 2 (h k+2 + \fnK h k ) 
and then from Proposition |4.1[ we get that 

+ -t Y -)c 2 (h k+2 + \fnK h k ) + ^-C x h k ~ x . 

VytihJ y/v 


\y-Vh\\z s 


<21 


(56) 


Similarly, 


|A — A*.||jsc' < ( 2 


: + — C 2 (h k+2 + Vn.Kh k ) + -^-C 1 h k ~ 1 . 
\fyoh ) Oh \/yo h 


From the last two estimates, we obtain the conclusion of the proposition. □ 

It remains now to deduce the convergence of the approximated solution yh for the L 2 (Qt) 
norm: this is done using the observability estimate ©■ Precisely, we write that (y — yh) solves 

L (y - yh) = -Ly h in Q T 
((y - yh), (y - Vh)t){Q) e h 
y-yh = 0 on e t . 

Therefore using <§. there exists a constant C(Cn,T,C 0 b s ) such that 

II V ~ 2/b|lL 2 (Q T ) ^ C(Cn iT) C obs )(\\y — yh\\ 2 L 2 ( qT ) + \\Ly h \\ 2 x) 
from which we deduce, in view of the definition of the norm Y, that 


\y — yh\\L 2 (Q T ) < C(Cn j T, C 0 b s ) max(l, )||y — yh\\z- 


(57) 


Eventually, by coupling (57) and Proposition 4.2 we obtain the following result : 


Theorem 4.1 (BFS element for N = 1 - Rate of convergence for the norm L 2 (Q t )) 
Assume that the hypothesis © holds. Let h > 0, let k G { 1,2} be a positive integer and let 77 < 1. 
Let (y, A) and (yh,^h) he the solution of & and (41\ ) respectively. If the solution (y , A) belongs 
to H k+2 (QT) x H u (Qt), then there exists two positives constant K = K(\\y\\ H k+ 2 ^Q T ^ || c llc 7 i(Q^T)’ 
IMI|.L°°(Q r )> Cq,t, C 0 bs), independent of h, such that 

\\y — Vh\\L 2 (Q T ) < /imax(l, ^) h -^(^ y /fi+ ^){h 3 + y/rjh) + l^J. (58) 
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Remark 11 Estimate (58) is not fully satisfactory as it depends on the constant Sh ■ In view of 


the complexity of both the constraint Ly = 0 and of the structure of the space Zh , the theoretical 
estimation of the constant Sh with respect to h is a difficult problem. However, as discussed at length 
in JTB , Section 2.1], Sh can be evaluated numerically for any h, as the solution of the following 
generalized eigenvalue problem (taking y = r, so that a r (y,y) is exactly ||y|||}: 


5 h = inf : B h A^ h Bl{\ h } = S J h {X h }, V{A h } £ R mh \ 


{ 0 } 


(59) 


where the matrix A ry h, Bh and Jh are defined in HD- 


Table [7] reports the values of Sh for r = 1 and r = h~ 2 for several values of h, T = 2, 
u> = (0.1, 0.3) and the BFS element. As in \12] where the boundary situation is considered with 
more details, these values suggests that, asymptotically in h, the constant S rt h behaves like : 


S r h ~ C r —.- as h —>■ 0 4 

Vr 

with C r > 0, a uniformly bounded constant w.r.t. h. 


(60) 


h 

7.01 x 10" 2 

3.53 x 10 -2 

1.76 x 10" 2 

8.83 x 10" 3 

r = 1 

3.58 

3.48 

3.42 

3.40 

(N 

1 

II 

2.53 x 10" 1 

1.23 x 10" 1 

6.05 x 10 -2 

3.01 x 10" 2 


Table 1: £ = 0: T = 2 - S Tt h for r = 1 and r = h 2 with respect to h. 


Consequently, in view of 60 the right hand side of the estimate (58) of \\y — yh\\L 2 (Q T ) behaves, 
taking r/ = r and r > 1 so that max(l, AC) = 1, like 

II V ~ VhWl^Qt) < Khf- 1 L/fh + 

and reaches its minimum for r = 1 /h, leading to \\y — yh\\L 2 (Q T ) A Kh k ~ 1 A. 


Eventually, when the space Zh is based on the HCT element, Theorem |4 .1 1 and Remark 11 still 
hold for k = 1. From 0 ch. VII, sec. 48, p. 295], we use that, for k S {0,1}, there exists a 
constant C 2 > 0 such 


\y- n z h ,r h {y)\\L 2 (Q T ) < c 2 h 


k +2 


\\y - B- Zh ,T h {y)\\H 2 (Q T ) <C 2 h k . 


(61) 


Then, we use that the error \\y — yh\\L 2 (Q T ) is again controlled by the error on the Lagrange 
multiplier A through the term d(X,A h ) in (|43| to conclude. 


4.2 Numerical approximation of the mixed formulation (19) 


We address the numerical approximation of the stabilized mixed formulation (19) with a £ (0,1) 
and r > 0. Let h be a real parameter. Let Zh and be two finite dimensional spaces such that 

Z h C Z, A h C A, V7i > 0. 


The problem (19) becomes : find ( yh,Xh ) £ Z h x A h solution of 


a r , a {y h ,y h ) + b a {X h ,y h ) = h, a (y h )> 

ba(Xh, Ph ) Ca(Xh, A^) ^2,a(A/j), 


^Vh ^ ^h 

VA h £ A h , 


(62) 


Proceeding as in the proof of [U Theorem 5.5.2], we first easily show that the following estimate 
holds . 
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Lemma 4.1 Let (y, A) £ Y x A be the solution of (19) and (yh,Xh) G Zh x fee tfce solution of 

,, Q || 2 , \ K \\ 2 


(62). Then we have, 




+ tt - iri t lly* -2/IU 


y h £Zh 

\b a \\ 2 a 2 d 2 \ ||T 

+ + TT _ m L IIAft. - A|| 


di 


6 , 


^ ' A^eA), 


(63) 


with ||o r , a || < max(l — a,p 1 r), ||6 Q || < max(7y l i 2 ,a). Parameters 9\ and 62 are defined in (20). 
Concerning the space Ah, since LXh should belong to L 2 (Qt), a natural choice is 

A h = {A e A(-, 0) = A t (-, 0) = 0}. (64) 


where Zh C Z is defined by (49). Then, using Lemma |4.l| and the estimate (55), we obtain the 
following result. 

Proposition 4.3 (BFS element for N = 1 - Rate of convergence - Stabilized formulation) 
Let h > 0, let k < 2 be a positive integer and let a £ (0,1). Let (y, A) and ( yh , A^) be the solution of 
(19) and (62) respectively. If (y, A) belongs to H k+ 2 (QT) x H k+ 2 (QT), then there exists a positive 
constant I\ = K(\\y\\ H k+ 2 ^Q T ^, ||c||( 7 i(QR), ||d|| L°°(Q T ): a UUl) independent of h, such that 


\\y — Vh\\z + || A — Aft || A < I\h k . 
In particular, arguing as in the previous section, we get 


(65) 


Theorem 4.2 (Rate of convergence for the norm L 2 (Q t ) Stabilized formulation) Assume 
that the hypothesis m holds. Let h > 0, let an integer k < 2. Let (y, A) and ( J/ft, Aft) be the so¬ 
lution of (19) and (62) respectively. If the solution (y, A) belongs to H k+ 2 {QT ) x H k+ 2 (QT), then 
there exist a positive constant K = K(\\y\\ H k+ 2 (Q T ), || M\h^(q t ), IMIcRQTp \\ d \\L°°(Q T ),a, r ,v) in¬ 
dependent of h such that 

h k 

\\V - Vh\\LRQ T ) < k—. ( 66 ) 


5 Numerical experiments 


We now report and discuss some numerical experiments corresponding to mixed formulation (41) 


and (62) for N = 1 and N = 2. 


5.1 One dimensional case (N = 1) 

We take f2 = (0,1). In order to check the convergence of the method, we consider explicit solutions 
of 0. We define the smooth initial condition (see my 

( yo(x) = 16x 2 (l — x) 2 , 

(EX1) \ WV , , , 16(0,1) 

1 Vi(x) = (3* - 4x 3 ) l ( o,o. 5 )(z) + ( 4 x - 12x 2 + 9x - 1) l (0 . 5il) (x), 

and / = 0. The corresponding solution of (JT|) with c = 1, d = 0 is given by 
2 /(x, t) = n ak cos(/c7rt) + -1— sin(fc7rt) J v / 2sin(fc7rx) 


fc> 0 


with 


Ofc = 


32i/2(7r 2 fc 2 — 12) fc 48\/2sin(7rfc/2) 


TT 5 k 5 


-((-l) fe -l), b k = 


7 r 4 fc 4 


, fc > 0. 
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We also define the initial data in Hq(LI) x L 2 (Ll) 

(EX2) y 0 (x) = 1 - \2x- 1|, yi(x) = 1 (i/ 3 , 2 / 3 ) 0*0. a; e (0,1) 
for which the Fourier coefficients are 

Ofc = —sin(7 rfc/2), 6*, = -^-(cos(7rfc/3) — cos(27rfc/3)), k > 0. 

7 r^k^ nk 

5.1.1 The cylindrical case: qx = to x (0 , T) 

We consider the case e = 0 described in Section[2] We take co = (0.1,0.3) and T = 2 for which the 
inequality © holds true. We consider the BFS finite element with uniform triangulation (each 
element K of the triangulation 7), is a rectangle of lengths Ax and At so that h = \J (Ax) 2 + (At) 2 ). 
We recall that the direct method amounts to solve, for any h, the linear system ( |48| . We use the 
LU decomposition method. Table [2] collects some norms with respect to h for the initial data 
(EX1) for r = 1 and for Ax = At. We observe a linear convergence for the variables yh, A h for 
the L 2 -norm: 


h yh\\Li(Q T) =o(/t i. 03)i Hy Vh\\m qT ) = 0 (fe o. 98); || A/i[ | L2{ ) = 0 (h 0 - 98 ). (67) 

II2/IU 2 (Qt) \\y\\L 2 (q T ) 

In agreement with Remark [lj since y Q b s is by construction the restriction to qx of a solution of ([Tj) , 
the sequence Ah, approximation of A, vanishes as h —> 0. The L 2 -norm of Lyh do also converges 
to 0 with h, with a lower rate: 

\\Ly h \\ L 2 (QT) =O(h 0n ). ( 68 ) 


h 

\\y-yh II i/i( q.j ) 

\\v\\lHq t ) 

\\v-yh\\ L 2 (qT) 

WvWl^Wt) 

\\ L Vh\\L 2 (Q T ) 

l|A/l||i2 ( Q T ) 

K 

card({A ft }) 
tt CG iterates 


7.01 x 10 " 2 
9.55 x 10 " 2 

8.35 x 10 " 2 

5.62 x 10 " 3 
2.67 x 10 " 5 
1.4 x 10 10 
861 
27 


3.53 x 10 ~ 2 
4.58 x 10 ~ 2 

4.28 x 10 " 2 

3.21 x 10 ~ 3 
1.37 x 10 ~ 5 
4.6 x 10 11 
3 321 
42 


1.76 x 10 " 2 
2.24 x 10 " 2 

2.16 x 10" 2 

1.78 x 10 " 3 
6.89 x 10 " 6 
1.3 x 10 13 
13 041 
70 


8.83 x 10 ~ 3 
1.10 x 10 -2 

1.09 x 10 ~ 2 

9.99 x 10 ^ 4 
3.44 x 10 ' 6 
4.2 x 10 14 
51 681 
96 


4.42 x 10 " 3 
5.52 x 10 " 3 

5.51 x 10 " 3 

8.54 x 10 " 4 
1.76 x 10 " 6 
1.3 x 10 16 
205 761 
90 


Table 2: Example EX1 - r= l- T = 2- ||j/||z, 2 ( gT ) = 5.95 x 10 2 - ||j/||l 2 (q t ) = 1-59 x 10 1 . 


We also check that the minimization of the functional J** introduced in Proposition 2.3 leads 
exactly to the same result: we recall that the minimization of the functional J** corresponds 
to the resolution of the associate mixed formulation by an iterative Uzawa type method. The 
minimization is done using a conjugate gradient algorithm ( we refer to jl2( Section 2.2] for the 
algorithm). Each iteration amounts to solve a linear system involving the matrix A r h which is 
sparse, symmetric and positive definite. The Cholesky method is used. The performance of the 
algorithm depends on the conditioning number of the operator P r : precisely, it is known that (see 
for instance ED. 


|A"-A|U 2(Qt) <2v / ^) 


1 

\J v(V r ) +1 


IIa 0 - A|| L 2(Q T ), Vn > 1 


where A minimizes J**. v(V r ) = H'Prll ||'P r x || denotes the condition number of the operator V r . As 
discussed in El Section 4.4], the conditioning number of V r restricted to A h C L 2 (Qx) behaves 
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asymptotically as C“ 2 /i _2 . Table[2]reports the number of iterations of the algorithm, initiated with 
A 0 = 0 in Qt■ We take e = 10“ 10 as a stopping threshold for the algorithm (the algorithm is stopped 
as soon as the norm of the residue g n given here by Ly n satisfies ||< 7 ™|| l 2 (q t ) < 10 in ||g°|| L 2(Q r )). 

Table [2] reports the number of iterates to reach convergence, with respect to h. We ob¬ 
serve that this number is sub-linear with respect to h, precisely, with respect to the dimension 
nih = card({Xh}) of the approximated problems. This renders this method very attractive from a 
numerical point of view. 

From Remark [6j we also check the convergence w.r.t. h when we assume from the beginning 
that the multiplier A vanishes (see Table [3]). This amounts to minimize the functional J r given by 
or, equivalently, to perform exactly one iteration of the conjugate gradient algorithm we have 
just discussed. With r = 1, we observe a weaker convergence : 

II y ~ yh\\L 2 {Q T ) _ 0^0.574) II y ~~ Uh\\L 2 (q T ) _ £,/^0.94\ 

\\y\ U 2 (q t ) ’ IMU 2 ojt) 

This example illustrates that the convergence of Lyu to 0 in the norm L 2 (0, T, iJ _1 (0,1)) is 
enough here to guarantee the convergence of the approximation y we get that h\\Lyh\\L 2 (Q T ) ~ 
ll^yh||L 2 (o,T;ff-i(o,i) = 0(h 03 ) while \\Ly h \\ L 2 (Q T) slightly increases. Obviously, in this specific 
situation, a larger r (acting as a penalty term) independent of h yields a lower \\Lyh\\L 2 (Q T ) norm. 


h 

7.01 x 10" 2 

3.53 x 10 -2 

1.76 x 10" 2 

8.83 x 10" 3 

4.42 x 10" 3 

\\y~yh IIe 2 (q t ) 

9.74 x 10" 2 

8.35 x 10 -2 

7.72 x 10" 3 

4.90 x 10" 2 

4.28 x 10" 2 

1.11 x 10" 2 

2.84 x 10" 2 

2.18 x 10" 2 

2.01 x 10" 2 

2.16 x 10" 2 

1.12 x 10" 2 

3.40 x 10" 2 

2.01 x 10" 2 

6.21 x 10" 3 

4.79 x 10" 2 

lly-y/.lli.2( 9T) 

H^Hl, 2 (q T ) 

\\ L yh\\L 2 (Q T ) 


Table 3: Example EX1 - r = l- T = 2- A fixed to zero. 


On the contrary, we check that the convergence to 0 of || y—Vh\\L 2 (Q T ) is lost when the inequality 
© is not satisfied: Table [4] collects the norms w.r.t. h for the same data except the value T = 1 
(for which the uniqueness of the solution is lost): we observe that ||y — yh\\L 2 (Q T ) increases as 
h —► 0. As an illustration of the loss of uniqueness, these value also yields to a larger conditioning 
number n of the matrix Arg¬ 


h 

7.01 x 10" 2 3.53 x 10" 2 1.76 x 10" 2 8.83 x 10" 3 4.42 x 10~ 3 

\\V-Vh\\ L 2 { Q T) 

1.21 x lO’ 1 1.08 x 10” 1 1.34 x 10" 1 2.42 x 10" 1 5.19 x lO” 1 

8.40 x 10" 2 4.34 x 10” 2 2.22 x 10" 2 1.12 x 10" 2 5.62 x 10~ 3 

5.62 x 10" 2 2.77 x 10” 2 2.63 x 10" 2 2.25 x 10" 2 2.15 x 10~ 2 

1.84 x 10" 5 9.48 x 10" 6 4.76 x 10" 6 2.38 x 10" 6 1.19 x 10” 6 

1.2 x 10 11 9.8 x 10 12 1.1 x 10 15 1.5 x 10 17 2.7 x 10 19 

\\v\\l 2 (q t ) 

\\v-Vh\\ L 2 WT) 

II«IIi, 2( 9t ) 

\\ L yh\\L 2 (Q T ) 

II^IU 2 (Qr) 

K 


uable 4: Example EX1 - r = l- T= l- \\y e x\\L 2 (q T ) = 4.21 x 10 2 - \\y e x\\L 2 (Q T ) = 1-12 x 10 A 

Similar conclusions hold with the less regular initial data (EX2). Numerical results are reported 
in Table [5] We still observe a linear convergence w.r.t. h of \\y — yh\[L 2 (Q T )-> 111/ — Dh\\L 2 (q T ) 
and || Xh\\L 2 (Q T )- One notable difference is that the convergence rate is weaker for the norm 
\\Lyh\\L 2 (Q T y- 

II Ly h \\ L2{QT) =O(h 0123 ). (70) 

Again, this is enough to guarantee the convergence of yh toward a solution of the wave equation: 
recall that then ||I/i/h||z, 2 ( 0)T; h-Rq,!)) = 0(h 1123 ). We also observe that the number of iterates 
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in the CG algorithm remains largely sub-linear but is slightly larger: precisely, we have j] iter 
= O(h~ on ). Table § illustrates the case T = 1 while Table [I] illustrates the minimization of J r 
(see [l7]), both for r = 1. 


h 

7.01 x 10" 2 3.53 x 10~ 2 1.76 x 10" 2 8.83 x 10” 3 4.42 x 10" 3 

\\y~Vh II L 2 (Qt ) 

\\v\\l 2 (q t ) 

\\y-yh\\ L 2( qT) 

\\y\\ L -2( qT ) 

\\ L yh\\L2(Q T ) 

II'Vi||l 2 (Qt) 

jj CG iterates 

1.01 x 10" 1 4.81 x 10~ 2 2.34 x 10" 2 1.15 x 10~ 2 5.68 x 10" 3 

1.34 x 10" 1 5.05 x KT 2 2.37 x 10" 2 1.16 x 10~ 2 5.80 x 10" 3 

7.18 x 10" 2 6.59 x KT 2 6.11 x 10" 2 5.55 x 10" 2 5.10 x 10" 2 

1.07 x 10" 4 4.70 x 10" 5 2.32 x 10~ 5 1.15 x 10“ 5 5.76 x 10" 6 

29 46 83 133 201 


Table 5: Example EX2 - r = 1 - T = 2 - \\y\\ L *( qT ) = 1.56 x 10" 1 - \\y\\ L RQ T ) = 4.14 x 1CT 1 . 


h 

7.01 x 10" 2 3.53 x 10" 2 1.76 x 10" 2 8.83 x 10" 3 4.42 x KT 3 

\\y~yh\\ l 2 (q t ) 

2.74 xlO" 1 4.15 x 10" 1 6.30 x 10" 1 1.21 2.62 

1.37 x 10’ 1 5.76 x 10” 2 2.89 x 10" 2 2.41 x 10" 2 7.76 x 10~ 3 

5.97 x 10" 2 4.96 x 10" 2 4.96 x 10" 2 4.52 x 10" 2 4.21 x KT 2 

4.97 x 10 -5 2.32 x 10" 5 1.15 x 10" 5 5.76 x 10" 5 2.87 x KT 6 

II y II l 2 (q-j- ) 
lly-y/.lli.2( 9Tl 
\\y\\L 2 (q T ) 

1 \Lyh\\L°-(Q T ) 
'MU 2 (Qt) 


Table 6: Example EX2 - r = 1 - T = 1 - \\y\\ L R qT ) = 1-104 x 1CT 1 - \\y\\ L RQ T ) = 2.93 x 10" 1 . 


h 

7.01 x 10" 2 

3.53 x 10 -2 

1.76 x 10" 2 

8.83 x 10 -3 

4.42 x 10" 3 

II y—yh\\ l 2 (q t ) 

1.02 x 10" 1 

1.34 x 10’ 1 

7.43 x 10" 2 

5.27 x 10" 2 

5.06 x 10" 2 

7.43 x 10" 2 

3.18 x 10" 2 

2.37 x 10" 2 

8.65 x 10 -2 

2.48 x 10" 2 

1.21 x 10" 2 

1.10 x 10" 1 

2.25 x 10" 2 

6.65 x 10" 3 

1.37 x 10" 2 

Ify II l 2 (q-j- ) 
lly-y/.lli.2( 9Tl 
\\y\\L 2 (q T ) 

\\Lyh | l 2 ( q t ) 


Table 7: Example EX2 - r = l- T = 2- A fixed to zero. 


We end this section with some numerical results for the stabilized mixed formulation (62). The 
main difference is that the multiplier A is approximated in a much richer space (see 641 leading 
to larger linear system. Table [8] consider the case of the example EX2 with T = 2 and a = 1/2. 
In order to compare with the formulation (41), we take again r = 1. We observe the convergence 
w.r.t h and obtain slightly better rates and constants than in Table [5] in particular, we have 
|| y — 2/^||L 2 (Q T )/||y||L 2 (Q T ) = O(h 1A0 ). This is partially due to the fact that the space A h used for 
the variable A h in ( [62]) is richer than the space A^ used in (41). However, for a = 0 leading to 
the non stabilized mixed formulation, the space Kh is too rich and produce poor result, while we 
obtain very similar results for any values of a in (0,1]. Finally, we also check that - in contrast 
with the mixed formulation (41) - the positive parameter r does not affect the numerical results. 

We also emphasize that this variational method which requires a finite element discretization of 
the time-space Qt is particularly well-adapted to mesh optimization. Still for the example EX2, 
Figure [I] depicts a sequence of four distinct meshes of Qt = (0,1) x (0,T): the sequence is initiated 
with a coarse and regularly distributed mesh. The three other meshes are successively obtained 
by local refinement based on the norm of the gradient of yh on each triangle of 7/. As expected, 
the refinement is concentrated around the lines of singularity of yh travelling in Qt, generated by 
the singularity of the initial position y 0 . The four meshes contain 792,2 108,7 902 and 14 717 
triangles respectively (see Table [9] . The results obtained using the reduced HCT finite element 
are reported in Table [9] 
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h 

7.01 x 10" 2 

3.53 x 10 -2 

1.76 x 10" 2 

8.83 x 10 -3 

4.42 x 10" 3 

\\V-Vh\\ L 2 { Q T) 

8.48 x 10" 2 

4.01 x 10" 2 

1.85 x 10" 2 

8.66 x 10" 3 

4.01 x 10" 3 

IMh 2 (Q T ) 






\\v-Vh\\ L 2 WT) 

2.80 x 10" 1 

7.26 x 10" 2 

2.61 x 10" 2 

1.12 x 10" 2 

5.05 x 10" 3 

ll»lll,2( 9r ) 






\\ L yh\\L2(Q T ) 

7.25 x 10~ 2 

6.59 x lO^ 2 

6.16 x 10 -2 

5.58 x 10~ 2 

5.08 x KT 2 

IIAft |i 2 (Q T ) 

4.11 x 10" 3 

2.04 x 10" 3 

1.49 x 10" 3 

1.01 x 10" 3 

7.37 x 10" 4 


Table 8: Example EX2 - r = 1 - T = 2 - a = 1/2 - ||2/||z, 2 (g- r ) = 5.95 x 10 2 - ||2 /||l 2 (q t ) = 
1.59 x 10" 1 . 



Figure 1: Iterative refinement of the triangular mesh over Qt with respect to the variable y. 
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Mesh number 

1 

2 

3 

4 

ft elements 

792 

2 108 

7 902 

14 717 

jt points 

429 

1 101 

4 041 

7 462 

\\y~yh ||^ 2 (q t ) 
IMI r 2 (Q r ) 

1.34 x 1CT 2 

8.69 x 10" 3 

6.01 x 10" 3 

5.9 x KT 3 

I^IU 2 (Qt) 

1.14 x 10" 5 

7.99 x 10 -6 

5.02 x 10" 6 

4.79 x 10" 6 


Table 9: Example (EX2) - Information concerning the meshes and approximation errors for mesh 
adaptation strategy. 


5.1.2 The non-cylindrical case 


We numerically illustrate the reconstruction of the state of the wave equation ([I]) from measure¬ 
ments y 0 bs which are available in domains qr C Qt non-constant in time (considered recently in 
[7] in a controllability context). Time dependent domains also appears for time under sampled 
observations (or measurements): we refer to [TTJ. In what follows we take T = 2 and qt to be one 
of the two following domains: 


:= < (, x , t) £ Qt such that 


x — 


3 1 
5T 


< — for every t £ (0, T) 


(71) 


<?T 




uG>! 


T T 

x * 4 * ’ ~2 




(72) 


These two pairs (T, q l T ) i = 1,2 satisfy the standard geometric optic condition: therefore, using 
[7], Proposition 2.1, inequality (j^]) holds true. Both domains q and q%, are displayed in Figure [ 2 ] 
with the coarsest of the meshes that are used for the numerical experiments in this section. 


2 - 


1.5- 


0.5 



0.5 

(a) 



Figure 2: Domain qp (a) and domain q^ (b) triangulated using some coarse meshes. 


We consider five levels of regular triangular meshes and use the reduced Hsieh-Clough-Tocher 
finite element. We illustrate our method on the reconstruction of the solution of the wave equation 


corresponding to initial data (EX2) considered in Section 5.1.1 
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Since domains q^ and q t satisfy the geometric optic condition, we obtain similar results as in 
the case qr = w x (0, T) studied in the previous section. More precisely, these results are reported 

for domain q ^ and respectively. 


in Table 10 and Table 


11 


h 

7.18 x 10' 

-2 

3.59 x KT 2 

1.79 

x 10 

-2 

9 x 10" 3 

4.5 x 10- 

3 

\\y~yh IIe 2 (q t ) 
H'ylli 2 (Q T ) 

2.02 x 10“ 

-2 

7.83 x KT 3 

3.32 

x 10 

-3 

1.36 x 10" 3 

6.27 x 10 

-4 

\\y-yh\\ L 2( qT) 

ll^lll, 2 (q T ) 

1.85 x 10“ 

-2 

6.69 x 10" 3 

2.40 

x 10 

-3 

1.03 x 10" 3 

4.56 x 10 

-4 

\\ L yh\\L 2 (Q T ) 

3.41 


3.78 

I 

1.15 


4.47 

4.76 


IIAft |l 2 (q t ) 

1.97 x 10“ 

-5 

7.03 x 10 -6 

1.70 

x 10 

-6 

4.14 x 10" 7 

1.10 x 10 

-7 

K 

1.18 x 10 

8 

1.84 x 10 9 

1.61 

x 10 

10 

1.75 x 10 11 

1.38 x 10 

12 

card({A/J) 

429 


1 633 

6 

369 


25 153 

99 969 


j) CG iterates 

108 


206 

392 


954 

2 009 



Table 10: Observation domain q)p. Example EX2 - r = 1 - T = 2- ||y||L 2 (q T ) = 2.75 x 10 1 - 
IMU 2 (Qt) = 5-87 x 10 _1 . 

Remark that the number of iterations needed for the conjugate gradient algorithm in order 
to achieve a residual smaller than 10“ 10 when we minimize the functional J** over is slightly 
larger than in the situations described in the previous section. 


h 

6.24 x 10" 2 

3.12 x KT 2 

1.56 x 10" 2 

7.8 x 10" 3 

3.9 x 10~ 3 

\\y~yh\\ l 2 (q t ) 
\\v\\l 2 (,q t ) 

\\y-Vh\\ L 2 {qT) 

\\y\\L->(q T ) 

\\ L yh\\L2(Q T ) 
Pbi | l 2 ( q t ) 

K 

card({A/J) 
tt CG iterates 

1.38 x 10" 2 

1.27 x 10" 2 

3.86 

6.37 x 10 -6 
2.02 x 10 s 
554 

141 

6.37 x 10" 3 

4.79 x 10" 3 

3.45 

1.65 x 10" 6 
2.62 x 10 9 

2 135 

331 

2.64 x 10" 3 

2.02 x 10" 3 

3.36 

3.88 x 10" 7 
2.05 x 10 10 

8 381 

720 

1.15 x 10" 3 

9.11 x 10" 4 

3.85 

9.74 x 10" 8 
1.61 x 10 11 
33 209 

1 446 

5.25 x 10" 4 

4.29 x 10" 4 

4.16 

2.90 x 10" 8 
1.32 x 10 12 
132 209 

3 318 


Table 11: Observation domain q^. Example EX2 - r= l- T = 2- ||y||L 2 (g T ) = 2.75 x 10 1 - 
\\y\\ LHQT) = 5.87 x IQ" 1 . 


The exact solution y corresponding to initial data (EX2) is displayed in Figure[3](a) using the 
third mesh of the domain in Figure [ 2 ] (b). Figure [3] (b) illustrates the solution yh of the mixed 
formulation (41), where the observation y 0 b s is obtained as the restriction of y to q^. 


5.2 Two-dimensional space case (N = 2) 

We now illustrate the method introduced in Section [2] in the two-dimensional case. The procedure 
is similar but a bit more involved on a computational point of view, since Qt is now a subset of 
R 3 . 

In order to approach the mixed-formulation |7]), we consider a mesh 7 ~h of the domain Qt = 
12 x (0,T) formed by triangular prisms. This mesh is obtained by extrapolating along the time 
axis a triangulation of the spatial domain 12. For an example in the case 12 = (0, l) 2 and T = 2 
see Figure [ 4 ] (b) and for an example in the case of non-rectangular domains 12 C R 2 see Figure [ 5 ] 
(b). For both examples, the extrapolation along the the time axis is uniform : the height of the 
prismatic elements At is constant. 
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(a) (b) 


Figure 3: Example (EX2) (a) Reference solution, (b) Solution reconstructed from the observation 
Hobs = y\q^- 



Figure 4: (a) Example of sets f 1 and u>. (b) Example of mesh for = (0, l) 2 and T = 2. 


Let Zh be the finite dimensional space defined as follows 

Z = i ' Ph = € c1 (Qt) ^\k x1K2 e ¥(K xy ),e\ Kt e Q(K t ) 1 

h \ Vh = o on £ t for every K = K XlX2 x K t G %■ J ’ 


(73) 


P(A';c l3;2 ) is the space of functions corresponding to the reduced Hsieh-Clough-Tocher (F1CT for 
short) C' 1 -element recalled in Section 4.1.1; Q(K t ) is a space of degree three polynomials on the 
interval K t of the form [tj, tj+ 1 ] defined uniquely by their value and the value of their first derivative 
at the point tj and tj+i ■ In other words, Y), is the finite element space obtained as a tensorial 
product between the reduced HCT finite element and cubic Hermite finite element. We check that 
on each element K = K XlX2 x K t , the function is determined uniquely in term of the values 
of T, k := {ip(ai),ip Xl (ai),<p X2 (ai),ipt(ai),(p Xu t(ai),<p X2tt (ai),i = 1,- • • ,6} at the six nodes a* of K. 
Therefore, dimSif = 36. 

Similarly, let be the finite dimensional space defined by 

ip h = ip(xi,x 2 )0{t) G C°(Q t ) ip\ KxiX2 G F 1 {K XlX2 ),d\ Kt G Qi(A't) 

<Ph = 0 on E t for every I< = K XlX2 x K t G 77, 















5 NUMERICAL EXPERIMENTS 


27 


where Pi (K XlX2 ) and Qi(K t ) are the spaces of degree one polynomials on the triangle K XlX2 and 
interval K t respectively. 

For any h, we check that Z k C Z and that C A. 


5.2.1 Wave equation in a square 

We first consider the case LI defined by the unit square and again some explicit solutions used in 
[8]. Precisely, we define the following smooth initial condition: 


or,') / yo(xi,x 2 ) = 256xlx?,{l- x 1 ) 2 (1-x 2 ) 2 , , 

(EX1_2D> 1 iftfe.x,) = (1 - |2x, - 1|)(1 - 2 x 2 - 1|) £ ° 

The corresponding solution of ([I]) with c = 1, d = 0 and / = 0 is given by : 

y(xi,X 2 ,t) = V' (a k i cos(nkit) + — sin (/q^t) ) sin(kTrx) sm(lny), 

where y k i = 7r \/k 2 + 1 2 for every k,l G Z* and 

i 0 (t r 2 fc 2 - 12 )(t t 2 1 2 - 12) t 


(75) 


(76) 


akl = 2 


bid = 


2 5 


7 T 10 k 5 l 5 
7T k . nl 


'((-!)* -1)((-1)‘ -1) 


7 Ck 2 P Sm 2 Sm 2 ' 


We also dehne the following initial data (yo,yi) G Hq(LI) x L' 2 (Ll): 

(EX2-2D) j = (1 - |2 *‘ “ ‘(I 1 - l 2 « ^ !D 

[ yi(xi,x 2 ) = i(|,| ) 2 ( :E i ) x v 

The Fourier coefficients of the corresponding solution are 

2 5 7rfc 7 tI 

° w = ^W sm T sin T 

l 


(£ 1 , 2 : 2 ) € fi. 


(77) 


frfcZ = 


r 2 /cZ 


( 7T k 

2nk\ 

( 7 tI 

2nl \ 

COS —« 

— cos- 

cos — 

— cos — 

V 3 

3 ) 

V 3 

3 ) 


In what follows, we consider w the subset of LI described in Figure [4] (a) and given by: 

w = ((0,0.2) x (0,1)) U ((0,1) x (0,0.2)). (78) 

It is easy to see that this choice of u> and T = 2 provide a domain q? = w x (0, T ) which satisfies 
the geometric optic condition, and, hence, inequality © holds. We consider 3 levels of meshes of 
Qt , labelled from 1 to 3 and containing the number of elements (prisms) and nodes listed in Table 

m 


Mesh Number 

1 

2 

3 

Number of elements 

5 320 

15 320 

42 230 

Number of nodes 

3 234 

8 799 

23 370 

A t 

0.2 

0.1 

0.05 


Table 12: Characteristics of the meshes used for Qt = (0, l) 2 x (0,2). 

For each of these meshes we solve the mixed formulation 0 with the term y 0 i, s appearing in 
the right-hand side obtained as the restriction to qr of the solution computed by (76) for initial 
data EX1—2D and EX2—2D. 
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Mesh number 

1 

2 

3 

\\y-vh\\ L 2 ( Q T) 

\\vI\l 2 (Q t ) 

4.58 x 10" 2 

3.18 x 10" 2 

1.38 x 10" 2 

\\Lyh\\LRQ T ) 

1.44 

1.05 

1.05 

II'MUaqt) 

2.87 x 10" 5 

1.36 x 10" 5 

7.34 x 10" 6 

ft CG iterates 

121 

180 

168 


Table 13: e = 0: Example EX1 2D - r = 1. 


Table [13| concerns the example EX1—2D. In this table we list the norm of the relative error 
between the exact solution y given by (761 and the solution y h of the mixed formulation 0, the 
L 2 norm of Lyh and the L 2 norm of the Lagrange multiplier A h- 

As theoretically stated in Remark [l] and observed in numerical experiments in the case N = 1 
(see, for instance, Table 13), the Lagrange multiplier A h vanishes as h —> 0. In Table 14 we display 
the results obtained by numerically solving the variational problem 0 obtained from the mixed 
formulation when A h = 0. 


Mesh number 

1 2 3 

\\y-vh\\ L 2 (QT) 
\\y\\L 2 (Q T ) 
\\ L Vh\\ l 2 (Q t ) 

7.05 x 10" 2 4.44 x 10~ 2 2.37 x 10" 2 

1.31 0.97 0.97 


Table 14: Example EX1 2D r = 1 - A fixed to zero. 


Tables |l5| and [l6| display the results obtained for the initial data specified by EX2 2D, for the 
solutions (yh, A h) of the mixed formulation and for the variational problem obtained when Xh = 0 
respectively. 


Mesh number 

1 

2 

3 

\\y-vh\\ L 2 (QT) 

\\y\\L 2 (Q T ) 

4.74 x 10" 2 

3.72 x 10" 2 

2.09 x 10" 2 

\\ L Vh\\ l 2 (Qt) 

1.18 

0.89 

1.06 

U 2 (Qt) 

3.21 x 10 -5 

1.46 x 10" 5 

1.17 x 10" 5 

jj CG iterates 

128 

191 

168 


Table 15: Example EX2 2D - r = 1. 

The results are similar for both examples. In both cases we observe a linear convergence of 
yh to y in the norm L 2 over Qt when h goes to zero. Similarly, the norm ||A/j||£ 2 (q t ) linearly 
decreases as h goes to zero. 


Mesh number 

1 2 3 

\\y-yh\\ L 2 (QT) 

l|-^2/h||L 2 (Q T ) 

6.75 x 10" 2 4.93 x 10~ 2 3.37 x 10" 2 

1.07 0.82 0.97 


Table 16: Example EX2—2D r = 1 - A fixed to zero. 


5.2.2 Wave equation in a non-rectangular domain of R 2 

Let C R 2 be a domain with a regular boundary and w a non-empty subset with regular boundary. 
An example of such a configuration is illustrated in Figure 0] (a). As in the previous section, we 
take T = 2 and we build a mesh formed by triangular prisms of the domain Qt = LI x (0, T). An 
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example of such a mesh associated to the domain L1 is displayed in Figure [5] (b). This mesh is 
composed by 17 934 nodes distributed in 32 140 prismatic elements (this mesh corresponds to the 
mesh number 2 described in Table 17). 



Figure 5: (a) Example of sets LI and u. (b) Example of mesh of the domain Qt- 


We consider three levels of meshes of the domain Qt formed by the number of prisms and 
containing the number of nodes reported in Table |17[ 


Mesh number 

1 

2 

3 

Number of elements 

5 730 

32 1400 

130 280 

Number of nodes 

3 432 

17 934 

69 864 

Height of elements (At) 

0.2 

0.1 

0.05 


Table 17: Characteristics of the three meshes associated with Qt- 


Comparing to the situation described in Subsection |5.2.1[ the eigenfunctions and eigenvectors 
of the Dirichlet Laplace operator defined on LI are not explicitly available here. Consequently, 
from a given set of initial data, we numerically solve the wave equation |l]) using a standard time¬ 
marching method, from which we can extract an observation on qt- Precisely, we use a P\ finite 
elements method in space coupled with a Newmark unconditionally stable scheme for the time 
discretization. Hence, we solve the wave equation on the same mesh which was extrapolated in 
time in order to obtain the mesh number 2 of Qt- This two-dimensional mesh contains 1 704 
nodes and 3 257 triangles. The time discretization step is At = 10~ 2 . We denote y h the solution 
obtained in this way for the initial data (yoiVi) L Hq(LI) x L 2 (Ll) given by 


f -Ay 0 = 10, in LI 
l 2/o=0, on dLl, 


(79) 


From y h we generate the observation y 0 b s as the restriction of y h to qr- Finally, from this 
observation we reconstruct yh as the solution of the mixed formulation (29) on each of the three 
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meshes 

meshes 


described in Table 17 Table [18] display some norms of yh and A h obtained for the three 
and illustrates again the convergence of the method. 


Mesh number 

1 

2 

3 

\\Vh~ Vh \\l 2 (Q t ) 

II Vh Hl 2 (Q t ) 

1.88 x HT 1 

8.04 x 10" 2 

7.11 x 10" 2 

\\Lyh\\L^(Q T ) 

3.21 

2.01 

1.57 

IIAft, |l 2(Q t ) 

8.26 x 10" 5 

3.62 x 10" 5 

2.84 x 10" 5 

jj CG iterates 

52 

167 

400 


Table 18: Initial data (yo,yi) given by 


r = 1. 


Figure[6](a) displays the solution y 0 of (79) and Figure[6](b) displays the initial position y h (•, 0) 
corresponding to the solution of our inverse problem. The error between these two functions is 
given by ||j/o — Uh(',0)\\ L 2/^ = 2.05 x 10 -2 which is consistent with the results reported in Table 



Figure 6: (a) Initial data yo given by ([79|. (b) Reconstructed initial data yh(-, 0). 


6 Concluding remarks and perspectives 

The mixed formulations we have introduced here in order to address inverse problems for the wave 
equation seems original. These formulations are nothing else than the Euler systems associated 
to least-squares type functionals and depend on both the state to be reconstruct and a Lagrange 
multiplier. This Lagrange multiplier is introduced to take into account the state constraint Ly—f = 
0 and turns out to be the controlled solution of a wave equation with the source term (y—y 0 bs) 1 q T - 
This approach, recently used in a controllability context in |12j . leads to a variational problem 
defined over time-space functional Hilbert spaces, without distinction between the time and the 
space variable. The main ingredient allowing to prove the well-posedness of the mixed formulation 
and therefore the reconstruction of the solution, is a generalized observability inequality, assuming 
here some geometric conditions on the observation zone. 

At the practical level, the discrete mixed time-space formulation is solved in a systematic way 
in the framework of the finite element theory. The approximation is conformal allowing to obtain 
the strong convergence of the approximation as the discretization parameters tends to zero. In 
particular, we emphasize that there is no need, contrarily to the classical approach, to prove some 










REFERENCES 


31 


uniform discrete observability inequality: we simply use the observability equality on the finite 
dimensional discrete space. The resolution amounts to solve a sparse symmetric linear system : 
the corresponding matrix can be preconditioned if necessary, and may be computed once for all 
as it does not depend on the observation y 0 b s - Eventually, the space-time discretization of the 
domain allows an adaptation of the mesh so as to reduce the computational cost and capture the 
main features of the solutions. Similarly, this space-time formulation is very appropriate to the 
non-cylindrical situation. 

In agreement with the theoretical convergence, the numerical experiments reported here display 
a very good behavior and robustness of the approach: the reconstructed approximate solution 
converges strongly to the solution of the wave equation associated to the available observation. 
Remark that from the continuous dependence of the solution with respect to the observation, the 
method is robust with respect to the possible noise on the data. 

As mentioned at the end of Section [3j additional assumption on the source term allows to 
determine uniquely the pair (y, /) from a partial measurement on qx or on a part Ey sufficiently 
large of the boundary. For instance, from [231 Theorem 2.1], assuming that the source term takes 
the form f(x,t ) = a(t)y(x) with er £ C 1 ([0, T]), er(0) / 0 and y £ H~ 1 ( fl), then the following 
holds: there exists a positive constant C such that 


INI 2 h-i 


(fi) < c 


dy 

dv 


L 2 (£t) 


+ 11 Ly-<j{t)n(x) 


I i 2 (Qi 


V(y, y) e S 


(80) 


where y solves |T]) with (yo,yi) = 0, c = 1 and ( Ex,T,Qx ) satisfies a geometric condition and S 
denotes an appropriate functional space. Using this inequality (similar to HD- we can study the 
mixed formulation associated to the Lagrangian from S x L 2 {Qx) R. defined by 


to fully reconstruct y and y from y Q b s and a. 

Eventually, since the mixed formulations rely essentially on a generalized observability inequal¬ 
ity, it may be employed to any other observable systems for which such property is available : we 
mention notably the parabolic case usually - in view of regularization property - badly conditioned 
and for which direct and robust methods are certainly very advantageous. We refer to |23| where 
this issue is investigated. 


dy 

dv 


dobs 


rJCSrr, i 


A (Ly — ay) dx dt 
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